"""
精确对角化，有上下自旋
"""

using JLD2
using Dates
#using LinearAlgebra

"""
统计一共有几个可能的状态
"""
function count_state_num(npart, nsite)
    num = 1
    for idx in 0:1:npart-1
        num = num * (nsite - idx)
    end
    den = 1
    for idx in 2:1:npart
        den = den * idx
    end
    return Int64((num // den)^2)
end


#println(count_state_num(2, 12))


"""
下一个状态
"""
function next_state(sta, maxint)
    sint = parse(Int64, sta, base=2)
    npart  = count_ones(sint)
    nsite = length(sta)
    if sint == maxint
        return lpad(repeat("1", npart), nsite, "0"), true
    end
    #println(npart)
    sint += 1
    while sint <= maxint
        if count_ones(sint) == npart
            break
        end
        sint += 1
    end
    return lpad(string(sint, base=2), nsite, "0"), false
end


#println(next_state("0011", parse(Int64, "1100", base=2)))
#println(next_state("1100", parse(Int64, "1100", base=2)))


"""
获得数组
"""
function get_state_array(npart, nsite)
    upsta = lpad(repeat("1", npart), nsite, "0")
    dnsta = lpad(repeat("1", npart), nsite, "0")
    res = Vector{String}(undef, count_state_num(npart, nsite))
    resinv = Dict{String, Int}()
    maxint = parse(Int64, rpad(repeat("1", npart), nsite, "0"), base=2)
    #println("开始状态 ", upsta * dnsta)
    for idx in 1:1:length(res)
        res[idx] = upsta * dnsta
        resinv[upsta*dnsta] = idx
        dnsta, flag = next_state(dnsta, maxint)
        #如果flag,dnsta就转回来了，也需要给upsta进位
        if flag
            upsta, flag = next_state(upsta, maxint)
        end
    end
    #println("开始状态 ", upsta * dnsta)
    return res, resinv
end


function have_bond(sta1, sta2, sic, sia)
    sig1, statmp = apply_annihilation(sta1, sia)
    if ismissing(statmp)
        return 0
    end
    sig2, statmp = apply_creation(statmp, sic)
    if ismissing(statmp) || statmp != sta2
        return 0
    end
    return sig1*sig2
end

"""
检测从第一个状态到第二个状态，有几个bond
"""
macro have_bond(sta1, sta2, bdc, bda)
    return esc(quote
        bnum = 0
        sta3 = $sta1[1:$bda-1] * "0" * $sta1[$bda+1:end]
        sta3 = sta3[1:$bdc-1] * "1" * sta3[$bdc+1:end]
        if sta3 != $sta2 || $sta1 == $sta2
            bnum = 0
        else
            bdl = min($bda, $bdc)
            bdr = max($bda, $bdc)
            sigp = 0
            #println($sta1, " ", bdl," ", bdr, " ", $sta1[bdl:bdr-1])
            for idx in bdl:1:bdr-1
                if $sta1[idx] == '1'
                    sigp += 1
                end
            end
            if $bda < $bdc
                sigp -= 1
            end
            bnum = (-1)^sigp
        end
        bnum
    end)
end


#get_state_array(2, 12)

#function test_bond()
#    println(@have_bond "0011" "0110" 2 3)
#    println(@have_bond "0011" "0110" 2 4)
#end

#test_bond()


function full_hamiltonian(npart, nsite, allstas, invstas, bonds, U)
    stacnt = length(allstas)
    ham = zeros(Float64, stacnt, stacnt)
    for coli in Base.OneTo(stacnt)
        colsta = allstas[coli]
        upsta, dnsta = colsta[1:nsite], colsta[nsite+1:end]
        # up hoppings
        for bnd in bonds
            sic, sia = bnd[1], bnd[2]
            sigup, upsta2 = apply_hoppings(upsta, sic, sia)
            if !ismissing(upsta2)
                ham[invstas[upsta2*dnsta], coli] = sigup*bnd[3]
            end
            #
            sic, sia = bnd[2], bnd[1]
            sigup, upsta2 = apply_hoppings(upsta, sic, sia)
            if !ismissing(upsta2)
                ham[invstas[upsta2*dnsta], coli] = sigup*bnd[3]
            end
        end
        # dn bonds
        for bnd in bonds
            sic, sia = bnd[1], bnd[2]
            sigdn, dnsta2 = apply_hoppings(dnsta, sic, sia)
            if !ismissing(dnsta2)
                ham[invstas[upsta*dnsta2], coli] = sigdn*bnd[3]
            end
            #
            sic, sia = bnd[2], bnd[1]
            sigdn, dnsta2 = apply_hoppings(dnsta, sic, sia)
            if !ismissing(dnsta2)
                ham[invstas[upsta*dnsta2], coli] = sigdn*bnd[3]
            end
        end

        # U项
        for st in Base.OneTo(nsite)
            if upsta[st] == '1' && dnsta[st] == '1'
                ham[coli, coli] += U
            end
        end
    end
    return ham
end


"""创建稀疏矩阵"""
function full_hmltCSC(npart, nsite, allstas, invstas, bonds, U)
    stacnt = length(allstas)
    cols = []
    rows = []
    vals::Vector{Float64} = []
    prtcnt = 1
    @info "start"*string(Dates.now())
    for coli in Base.OneTo(stacnt)
        if coli == 10^prtcnt
            @info string(10^prtcnt)*string(Dates.now())
            prtcnt += 1
        end
        colsta = allstas[coli]
        upsta, dnsta = colsta[1:nsite], colsta[nsite+1:end]
        # up hoppings
        for bnd in bonds
            sic, sia = bnd[1], bnd[2]
            sigup, upsta2 = apply_hoppings(upsta, sic, sia)
            if !ismissing(upsta2)
                push!(rows, invstas[upsta2*dnsta])
                push!(cols, coli)
                push!(vals, sigup*bnd[3])
                #ham[invstas[upsta2*dnsta], coli] = sigup*bnd[3]
            end
            #
            sic, sia = bnd[2], bnd[1]
            sigup, upsta2 = apply_hoppings(upsta, sic, sia)
            if !ismissing(upsta2)
                push!(rows, invstas[upsta2*dnsta])
                push!(cols, coli)
                push!(vals, sigup*bnd[3])
                #ham[invstas[upsta2*dnsta], coli] = sigup*bnd[3]
            end
        end
        # dn bonds
        for bnd in bonds
            sic, sia = bnd[1], bnd[2]
            sigdn, dnsta2 = apply_hoppings(dnsta, sic, sia)
            if !ismissing(dnsta2)
                push!(rows, invstas[upsta*dnsta2])
                push!(cols, coli)
                push!(vals, sigdn*bnd[3])
                #ham[invstas[upsta*dnsta2], coli] = sigdn*bnd[3]
            end
            #
            sic, sia = bnd[2], bnd[1]
            sigdn, dnsta2 = apply_hoppings(dnsta, sic, sia)
            if !ismissing(dnsta2)
                push!(rows, invstas[upsta*dnsta2])
                push!(cols, coli)
                push!(vals, sigdn*bnd[3])
                #ham[invstas[upsta*dnsta2], coli] = sigdn*bnd[3]
            end
        end

        # U项
        for st in Base.OneTo(nsite)
            valtmp = 0.0
            if upsta[st] == '1' && dnsta[st] == '1'
                valtmp += U
                #ham[coli, coli] += U
            end
            push!(rows, coli)
            push!(cols, coli)
            push!(vals, valtmp)
        end
    end
    return sparse(rows, cols, vals)
end


"""
获取哈密顿量
"""
function get_ham(npart, nsite, bonds, U)
    staarr = get_state_array(npart, nsite)
    #println(staarr)
    ham = zeros(ComplexF32, length(staarr), length(staarr))
    for idx2 in 1:1:length(staarr)
    for idx1 in 1:1:length(staarr)
        sta1 = staarr[idx1]
        sta2 = staarr[idx2]
        upsta1, dnsta1 = sta1[1:nsite], sta1[nsite+1:end]
        upsta2, dnsta2 = sta2[1:nsite], sta2[nsite+1:end]
        val = complex(0., 0.)
        #println(sta1, " ", sta2, " ",@have_bond upsta1 upsta2 bonds[1][1] bonds[1][2])
        for bnd in bonds
            #后面的bnd是消灭
            #up的部分
            if dnsta1 == dnsta2
                coef = @have_bond upsta2 upsta1 bnd[1] bnd[2]
                val += coef * bnd[3]
            end
            #dn的部分
            if upsta1 == upsta2
                coef = @have_bond dnsta2 dnsta1 bnd[1] bnd[2]
                val += coef * bnd[3]
            end
        end
        ham[idx1, idx2] += val
        ham[idx2, idx1] += conj(val) 
    end
    end

    #onsite U
    for idx in 1:1:length(staarr)
        sta = staarr[idx]
        upsta, dnsta = sta[1:nsite], sta[nsite+1:end]
        for st in 1:1:nsite
            if upsta[st] == '1' && dnsta[st] == '1'
                ham[idx, idx] += U
            end
        end
    end

    #V
    #for idx in 1:1:length(staarr)
    #    sta = staarr[idx]
    #    upsta, dnsta = sta[1:nsite], sta[nsite+1:end]
    #    #前面一半是第一个轨道
    #    for st1 in 1:1:Int64(nsite//2)
    #        n1 = 0.
    #        if upsta[st1] == '1'
    #            n1 += 1
    #        end
    #        if dnsta[st1] == '1'
    #            n1 += 1
    #        end
    #        st2 = st1 + Int64(nsite//2)
    #        n2 = 0.
    #        if upsta[st2] == '1'
    #            n2 += 1
    #        end
    #        if dnsta[st2] == '1'
    #            n2 += 1
    #        end
    #        ham[idx, idx] += 0.0 * n1 * n2
    #    end
    #end
    return ham
end


function test_ham1()
    println("a")
    println(@have_bond "100" "010" 2 1)
    bnds = [(1, 2, -1), (2, 1, -1)]#, (2, 4, -1), (1, 3, -1), (3, 4, -1)]

    ham = get_ham(1, 3, bnds, 8.0)

    sigham =  zeros(ComplexF32, 3, 3)
    for b in bnds
        sigham[b[1], b[2]] += b[3]
        #sigham[b[2], b[1]] += b[3]
    end
    println(real(sigham))
    println(eigvals(sigham))

    println(real(ham))
    evals, evecs = eigen(ham)

    println(evals)
    println(dot(evecs[:, 1], evecs[:, 1]))
end

#test_ham1()

#ham = get_ham(2, 12, [
#    (1, 2, -1),
#    (1, 4, -1),  
#    (1, 6, -1),
#    (2, 3, -1),
#    (2, 5, -1),
#    (3, 4, -1),
#    (3, 6, -1),
#    (4, 5, -1),
#    (5, 6, -1),
#    (7,  8, -1),
#    (7, 10, -1),
#    (7, 12, -1),
#    (8,  9, -1),
#    (8, 11, -1),
#    (9, 10, -1),
#    (9, 12, -1),
#    (10, 11, -1),
#    (11, 12, -1),
#    (1, 3, -complex(0.50,-0.50)), 
#    (1, 5, -complex(0.50,0.50 )),
#    (2, 4, -complex(0.50,0.50 )),
#    (2, 6, -complex(0.50,-0.50)), 
#    (3, 5, -complex(0.50,-0.50)), 
#    (4, 6, -complex(0.50,0.50)),  
#    (7, 9, -complex(0.50,0.50)),  
#    (7,  11, -complex(0.50, -0.50)),
#    (8,  10, -complex(0.50, -0.50 )),
#    (8,  12, -complex(0.50, 0.50 )),
#    (9,  11, -complex(0.50, 0.50 )),
#    (10,  12, -complex(0.50, -0.50)),
#], 8.0)

#println(eigvals(ham))


"""
作用hopping项
"""
function apply_hoppings(sta, bdc, bda)
    #消灭一个
    if sta[bda] == '0'
        return 0, missing
    end
    sigpow = 0
    for sts in sta[1:bda-1]
        if sts == '1'
            sigpow += 1
        end
    end
    statmp = sta[1:bda-1]*"0"*sta[bda+1:end]
    #产生一个
    if statmp[bdc] == '1'
        return 0, missing
    end
    for sts in statmp[1:bdc-1]
        if sts == '1'
            sigpow += 1
        end
    end
    return (-1)^sigpow, statmp[1:bdc-1]*"1"*statmp[bdc+1:end]
end


"""
作用哈密顿量
"""
function apply_hmlt!(u_tmp, ν_im1, nsite, staarr, stadic, bonds, U)
    u_tmp .= 0
    for idx in Base.OneTo(length(staarr))
        sta = staarr[idx]
        upsta = sta[1:nsite]
        dnsta = sta[nsite+1:end]
        #作用hopping项（上）
        for bnd in bonds
            sgn, sta2 = apply_hoppings(upsta, bnd[1], bnd[2])
            tpara = bnd[3]
            if !ismissing(sta2)
                #println(stadic[hres[2]*dnsta])
                u_tmp[stadic[sta2*dnsta]] += tpara * sgn * ν_im1[idx]
            end
            sgn, sta2 = apply_hoppings(upsta, bnd[2], bnd[1])
            #println(upsta, " ", bnd[1], " ", bnd[2], " ", hres)
            if !ismissing(sta2)
                u_tmp[stadic[sta2*dnsta]] += tpara * sgn * ν_im1[idx]
            end
        end
        #作用hopping项（下）
        for bnd in bonds
            sgn, sta2 = apply_hoppings(dnsta, bnd[1], bnd[2])
            tpara = bnd[3]
            if !ismissing(sta2)
                u_tmp[stadic[upsta*sta2]] += tpara * sgn * ν_im1[idx]
            end
            sgn, sta2 = apply_hoppings(dnsta, bnd[2], bnd[1])
            if !ismissing(sta2)
                u_tmp[stadic[upsta*sta2]] += tpara * sgn * ν_im1[idx]
            end
        end
        #
        # U
        for st in Base.OneTo(nsite)
            if upsta[st] == '1' && dnsta[st] == '1'
                u_tmp[idx] += U * ν_im1[idx]
            end
        end
    end
end


"""
构造lanzcos矩阵
"""
function lanzcos_iteration(npart, nsite, staarr, stadic, bonds, U)
    #
    if isfile("nutmp.jld")
        rm("nutmp.jld")
    end
    handle = jldopen("nutmp.jld", "w")
    #
    #staarr, stadic = get_state_array(npart, nsite)
    #println(stadic)
    #
    ν_im1 = Vector{ComplexF64}(undef, length(staarr))
    ν_i = Vector{ComplexF64}(undef, length(staarr))
    ν_ip1 = Vector{ComplexF64}(undef, length(staarr))
    #
    u_tmp = Vector{ComplexF64}(undef, length(staarr))#H*ν_im1
    #
    α_0 = Float64
    α_arr = Float64[]
    β_arr = Float64[]
    #初始的波函数ν_0, 归一化的
    ν_im1 .= rand(ComplexF64, length(staarr))
    ν_im1 = ν_im1 / real(sqrt(dot(ν_im1, ν_im1)))
    write(handle, "0", ν_im1)
    #通过ν_0构造ν_1
    apply_hmlt!(u_tmp, ν_im1, nsite, staarr, stadic, bonds, U)    
    #println(u_tmp)
    α_0 = real(dot(ν_im1, u_tmp))
    push!(α_arr, α_0)
    #println(α_0)
    #正交化
    u_tmp .= u_tmp - α_0 * ν_im1
    #归一化
    β = real(sqrt(dot(u_tmp, u_tmp)))
    push!(β_arr, β)
    ν_i .= u_tmp / β
    write(handle, "1", ν_i)
    #上面计算完了ν_1， β_1
    apply_hmlt!(u_tmp, ν_i, nsite, staarr, stadic, bonds, U)
    α = real(dot(ν_i, u_tmp))
    push!(α_arr, α)
    #println("ini ", dot(ν_i, ν_im1))
    #
    tmp_eng = α_0
    #开始迭代
    idx = 1
    while true
        #之前计算alpha的时候计算过了
        #apply_hmlt!(u_tmp, ν_i, nsite, staarr, stadic, bonds)
        u_tmp .= u_tmp .- β_arr[end]*ν_im1 .- α_arr[end]*ν_i
        β = real(sqrt(dot(u_tmp, u_tmp)))
        push!(β_arr, β)
        ν_ip1 .= u_tmp / β
        #计算 H ν_ip1
        apply_hmlt!(u_tmp, ν_ip1, nsite, staarr, stadic, bonds, U)
        α = real(dot(ν_ip1, u_tmp))
        push!(α_arr, α)
        #println(dot(ν_ip1, ν_im1), " ", dot(ν_ip1, ν_i))
        #向前推进一个
        ν_im1 .= ν_i
        ν_i .= ν_ip1
        idx += 1
        write(handle, string(idx), ν_i)
        #判断结果是否收敛
        if mod(idx, 5) == 0
            mat = SymTridiagonal(α_arr, β_arr)
            gndeng = eigmin(mat)
            if abs(tmp_eng - gndeng) < 1e-6
                tmp_eng = gndeng
                break
            end
            tmp_eng = gndeng
        end
        if length(α_arr) == length(staarr)
            mat = SymTridiagonal(α_arr, β_arr)
            gndeng = eigmin(mat)
            tmp_eng = gndeng
            break
        end
    end
    println("在第"*string(length(α_arr))*"处结束")
    #println(α_arr)
    #println(β_arr)
    #println(mat)
    #println(eigvals(mat))
    #ham = get_ham(npart, nsite, bonds, U)
    #println(eigvals(ham))
    #println(tmps[1], "aaa ", tmps[2])
    #println(dot(tmps[1], tmps[2]))
    #println(dot(tmps[1], tmps[end-1]))
    #println()
    #
    #println(ham)
    #println(ν_i)
    #println(ham * ν_im1)
    #保存基态
    println(tmp_eng)
    mat = SymTridiagonal(α_arr, β_arr)
    gndeng = tmp_eng#eigmin(mat)
    gndvec = eigvecs(mat, [tmp_eng])
    println(gndeng)
    #println(gndvec)
    write(handle, "gndeng", gndeng)
    write(handle, "gndvec", gndvec)
    close(handle)
end


"""
从结果重新构造
"""
function reconstruct_ground_state(npart, nsite, staarr, stadic, bonds, U)
    handle = jldopen("nutmp.jld", "r")
    rvec = read(handle, "gndvec")
    #println(rvec)
    #println(dot(rvec, rvec))
    #println(length(rvec))
    #println(rvec)
    resultvec = zeros(ComplexF64, length(staarr))
    chkvec = read(handle, "0")
    for idx = 1:1:length(rvec)
        tmpvec = read(handle, string(idx-1))
        #println(sum(tmpvec), " ", rvec[idx])
        #查看正交性
        println(dot(tmpvec, tmpvec), " ", dot(tmpvec, chkvec))
        resultvec .+= rvec[idx] * tmpvec
    end
    #println(resultvec)
    println(dot(resultvec, resultvec))
    tmpvec = Vector{ComplexF64}(undef, length(staarr))
    apply_hmlt!(tmpvec, resultvec, nsite, staarr, stadic, bonds, U)
    engcheck = dot(resultvec, tmpvec)
    gndeng = read(handle, "gndeng")
    println(engcheck)
    if abs(engcheck - gndeng) > 1e-5
        throw(error("low precision"))
    end
    close(handle)
    return resultvec
end


"""
计算自旋关联
"""
function spin_corr(npart, nsite, gndvec, st1, st2)
    staarr = get_state_array(npart, nsite)
    corr = 0.
    for idx = 1:1:length(staarr)
        sta = staarr[idx]
        upsta = sta[1:nsite]
        dnsta = sta[nsite+1:end]
        szst1 = 0.
        if upsta[st1] == '1'
            szst1 += 1
        end
        if dnsta[st1] == '1'
            szst1 -= 1
        end
        szst2 = 0.
        if upsta[st2] == '1'
            szst2 += 1
        end
        if dnsta[st2] == '1'
            szst2 -= 1
        end
        corr += adjoint(gndvec[idx]) * gndvec[idx] * szst1 * szst2
    end
    return corr
end


"""作用一个消灭算符"""
function apply_annihilation(sta, sti)
    if sta[sti] == '0'
        return 0, missing
    end
    sigpow = 0
    for sts in sta[1:sti-1]
        if sts == '1'
            sigpow += 1
        end
    end
    statmp = sta[1:sti-1]*"0"*sta[sti+1:end]
    return (-1)^sigpow, statmp
end


"""作用产生算符"""
function apply_creation(sta, sti)
    if sta[sti] == '1'
        return 0, missing
    end
    sigpow = 0
    for sts in sta[1:sti-1]
        if sts == '1'
            sigpow += 1
        end
    end
    statmp = sta[1:sti-1]*"1"*sta[sti+1:end]
    return (-1)^sigpow, statmp
end


"""
计算<v| 9 c^+_1,2 c^+_3,4 c_5,6 c_7,8 |v>
"""
function full_ppCSC(npart, nsite, allstas, invstas, ppbonds, U)
    stacnt = length(allstas)
    cols = []
    rows = []
    vals::Vector{ComplexF64} = []
    prtcnt = 1
    @info "start"*string(Dates.now())
    for coli in Base.OneTo(stacnt)
        if coli == 10^prtcnt
            @info string(10^prtcnt)*string(Dates.now())
            prtcnt += 1
        end
        colsta = allstas[coli]
        #每一个bnd会push一个 coli， 一个rowi和一个vals
        for bnd in ppbonds
            upsta, dnsta = colsta[1:nsite], colsta[nsite+1:end]
            tsig = 1.0
            #消灭第一个
            if bnd[8] == :up
                sig, upsta = apply_annihilation(upsta, bnd[7])
            elseif bnd[8] == :dn
                sig, dnsta = apply_annihilation(dnsta, bnd[7])
                sig *= (-1)^count_ones(parse(Int, upsta; base=2))
            else
                @assert false
            end
            if ismissing(upsta) || ismissing(dnsta)
                continue
            end
            tsig *= sig
            #消灭第二个
            if bnd[6] == :up
                sig, upsta = apply_annihilation(upsta, bnd[5])
            elseif bnd[6] == :dn
                sig, dnsta = apply_annihilation(dnsta, bnd[5])
                sig *= (-1)^count_ones(parse(Int, upsta; base=2))
            else
                @assert false
            end
            if ismissing(upsta) || ismissing(dnsta)
                continue
            end
            tsig *= sig
            #产生第一个
            if bnd[4] == :up
                sig, upsta = apply_creation(upsta, bnd[3])
            elseif bnd[4] == :dn
                sig, dnsta = apply_creation(dnsta, bnd[3])
                sig *= (-1)^count_ones(parse(Int, upsta; base=2))
            else
                @assert false
            end
            if ismissing(upsta) || ismissing(dnsta)
                continue
            end
            tsig *= sig
            #产生第二个
            if bnd[2] == :up
                sig, upsta = apply_creation(upsta, bnd[1])
            elseif bnd[2] == :dn
                sig, dnsta = apply_creation(dnsta, bnd[1])
                sig *= (-1)^count_ones(parse(Int, upsta; base=2))
            else
                @assert false
            end
            if ismissing(upsta) || ismissing(dnsta)
                continue
            end
            tsig *= sig
            push!(rows, invstas[upsta*dnsta])
            push!(cols, coli)
            push!(vals, tsig*bnd[9])
        end
    end
    push!(rows, stacnt)
    push!(cols, stacnt)
    push!(vals, 0.0)
    #println([(r,c,v) for (r,c,v) in zip(rows, cols, Val)])
    return sparse(rows, cols, vals)
end


"""
计算所有的关联函数
"""
function gfc_at_ij(i, j, sp, npart, nsite, allstas, invstas)
    stacnt = length(allstas)
    cols = []
    rows = []
    vals::Vector{ComplexF64} = []
    prtcnt = 1
    #@info "start"*string(Dates.now())
    for coli in Base.OneTo(stacnt)
        if coli == 10^prtcnt
            #@info string(10^prtcnt)*string(Dates.now())
            prtcnt += 1
        end
        colsta = allstas[coli]
        upsta, dnsta = colsta[1:nsite], colsta[nsite+1:end]
        if sp == :up
            sig, finsta = apply_hoppings(upsta, i, j)
            finsta = finsta * dnsta
        elseif sp == :dn
            sig, finsta = apply_hoppings(dnsta, i, j)
            finsta = upsta * finsta
        else
            @error "up or dn"
        end
        if !ismissing(finsta)
            push!(rows, invstas[finsta])
            push!(cols, coli)
            push!(vals, sig)
        end
    end
    push!(rows, stacnt)
    push!(cols, stacnt)
    push!(vals, 0.0)
    #println([(r,c,v) for (r,c,v) in zip(rows, cols, Val)])
    return sparse(rows, cols, vals)
end


"""
计算超导关联
"""
function spcd_corr(npart, nsite, gndvec, ord1, ord2)
    staarr = get_state_array(npart, nsite)
    stadic = Dict()
    for sta in staarr
        stadic[sta] = length(stadic) + 1
    end
    corr = 0.
    for idx = 1:1:length(staarr)
        sta = staarr[idx]
        #sta = "0000011000010010"
        upsta = sta[1:nsite]
        dnsta = sta[nsite+1:end]
        #println(sta)
        for bnd1 in ord1; for bnd2 in ord2
            #bnd1(ud) bnd2(du)
            #消灭两个
            ares1 = apply_annihilation(upsta, bnd2[1])
            if ismissing(ares1)
                ares1 = (0, upsta)
            end
            ares2 = apply_annihilation(dnsta, bnd2[2])
            if ismissing(ares2)
                ares2 = (0, dnsta)
            end
            #产生两个
            cres2 = apply_creation(ares2[2], bnd1[2])
            if ismissing(cres2)
                cres2 = (0, ares2[2])
            end
            cres1 = apply_creation(ares1[2], bnd1[1])
            if ismissing(cres1)
                cres1 = (0, ares1[2])
            end
            finalsta = cres1[2]*cres2[2]
            if finalsta in keys(stadic)
                finalvid = stadic[finalsta]
                #println("1 ", finalsta)
                #println(cres1[1], " ", cres2[1], " ", ares1[1], " ", ares2[1])
                corr += cres1[1] * cres2[1] * ares1[1] * ares2[1] *
                adjoint(gndvec[finalvid]) * gndvec[idx]
                #throw(error("a"))
            end
            #-bnd1(du) bnd2(du)
            #消灭两个
            ares1 = apply_annihilation(upsta, bnd2[1])
            if ismissing(ares1)
                ares1 = (0, upsta)
            end
            ares2 = apply_annihilation(dnsta, bnd2[2])
            if ismissing(ares2)
                ares2 = (0, dnsta)
            end
            #产生两个
            cres2 = apply_creation(ares1[2], bnd1[2])
            if ismissing(cres2)
                cres2 = (0, ares1[2])
            end
            cres1 = apply_creation(ares2[2], bnd1[1])
            if ismissing(cres1)
                cres1 = (0, ares2[2])
            end
            finalsta = cres2[2]*cres1[2]
            #println("2 ", finalsta)
            if finalsta in keys(stadic)
                finalvid = stadic[finalsta]
                #println("2 ", finalsta)
                corr -= cres1[1] * cres2[1] * ares1[1] * ares2[1] *
                adjoint(gndvec[finalvid]) * gndvec[idx]
            end
            #bnd1(du) bnd2(ud)
            #消灭两个
            ares1 = apply_annihilation(dnsta, bnd2[1])
            if ismissing(ares1)
                ares1 = (0, dnsta)
            end
            ares2 = apply_annihilation(upsta, bnd2[2])
            if ismissing(ares2)
                ares2 = (0, upsta)
            end
            #产生两个
            cres2 = apply_creation(ares2[2], bnd1[2])
            if ismissing(cres2)
                cres2 = (0, ares2[2])
            end
            cres1 = apply_creation(ares1[2], bnd1[1])
            if ismissing(cres1)
                cres1 = (0, ares1[2])
            end
            finalsta = cres2[2]*cres1[2]
            if finalsta in keys(stadic)
                finalvid = stadic[finalsta]
                #println("3 ", finalsta)
                corr += cres1[1] * cres2[1] * ares1[1] * ares2[1] *
                adjoint(gndvec[finalvid]) * gndvec[idx]
            end
            #-bnd1(ud) bnd2(ud)
            #消灭两个
            ares1 = apply_annihilation(dnsta, bnd2[1])
            if ismissing(ares1)
                ares1 = (0, dnsta)
            end
            ares2 = apply_annihilation(upsta, bnd2[2])
            if ismissing(ares2)
                ares2 = (0, upsta)
            end
            #产生两个
            cres2 = apply_creation(ares1[2], bnd1[2])
            if ismissing(cres2)
                cres2 = (0, ares1[2])
            end
            cres1 = apply_creation(ares2[2], bnd1[1])
            if ismissing(cres1)
                cres1 = (0, ares2[2])
            end
            finalsta = cres1[2]*cres2[2]
            if finalsta in keys(stadic)
                finalvid = stadic[finalsta]
                #println("4 ", finalsta)
                corr -= cres1[1] * cres2[1] * ares1[1] * ares2[1] *
                adjoint(gndvec[finalvid]) * gndvec[idx]
            end
        end; end
        #throw(error("a"))
    end
    return corr
end


#lanzcos_iteration(1, 6, 
#[(1, 2, -1.1), (1, 3, -0.9), (2, 3, -1.2), (3, 4, -0.95)], 
#1.23)

#staarr = get_state_array(2, 4)
#stadic = Dict()
#for sta in staarr
#    stadic[sta] = length(stadic) + 1
#end
#
#
#ham = get_ham(2, 4, [(1, 2, -1)], 0.0)
#ν_im1 = rand(ComplexF64, length(staarr))
#ν_im1 = ν_im1 / dot(ν_im1, ν_im1)
#println(ν_im1)
#
#println(ham*ν_im1)
#
#u_tmp = Vector{ComplexF64}(undef, length(staarr))
#apply_hmlt!(u_tmp, ν_im1, 4, staarr, stadic, [(1, 2, -1)], 0.0)
#println(u_tmp)

#function run()
#    println(@have_bond "1001" "0101" 1 2)
#    println(@have_bond "0101" "1001" 1 2)
#    ham = get_ham(1, 2, [(1, 2, -1)])
#    println(ham)
#end
#run()
#println(apply_hoppings("0101", 1, 4))

